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We present a general numerical scheme for the practical implementation of statistical moment clo- 
sures suitable for modeling complex, large-scale, nonlinear systems. Building on recently developed 
equation-free methods, this approach numerically integrates the closure dynamics, the equations of 
which may not even be available in closed form. Although closure dynamics introduce statistical 
assumptions of unknown validity, they can have significant computational advantages as they typi- 
cally have fewer degrees of freedom and may be much less stiff than the original detailed model. The 
closure method can in principle be applied to a wide class of nonlinear problems, including strongly- 
coupled systems (either deterministic or stochastic) for which there may be no scale separation. We 
demonstrate the equation-free approach for implementing entropy-based Eyink-Levermore closures 
on a nonlinear stochastic partial differential equation. 

PACS numbers: 



INTRODUCTION 

Accurate, fast simulations of complex, large-scale, non- 
linear systems remain a challenge for computational sci- 
ence and engineering, despite extraordinary advances in 
computing power. Examples range from molecular dy- 
namics simulations of proteins [l[ , Q and glasses 3 1 , to 
stochastic simulations of cellular biochemistry 0, 5|, to 
global-scale, geophysical fluid dynamics [|| . Often for the 
systems under consideration there is no obvious scale sep- 
aration, and their many degrees of freedom are strongly 
coupled. The complex and multiscale nature of these pro- 
cesses therefore makes them extremely difficult to model 
numerically. To make matters worse, one is often in- 
terested not in a single, time-dependent solution of the 
equations governing these processes, but rather in ensem- 
bles of solutions consisting of multiple realizations (e.g., 
sampling noise, initial conditions, and/or uncertain pa- 
rameters). Often real-time answers are needed (e.g., for 
control, tracking, filtering). These demands can easily 
exceed the computational resources available not only 
now but also for the foreseeable future. 

In principle, all statistical information for the problem 
under investigation is contained in solutions to the Liou- 
ville (if deterministic) /Kolmogorov (if stochastic) equa- 
tions. These are partial differential equations in a state 
space of high (possibly infinite) dimension. A straightfor- 
ward discretization of the Liouville / Kolmogorov equa- 
tions is therefore impractical. An ensemble approach to 



solving these equations can be taken; however, quite of- 
ten, the practical application of the ensemble approach 
is also problematic. Generating a sufficient number of 
independent samples for statistical convergence can be a 
challenge. For some problems, computing even one real- 
ization may be prohibitive. 

The traditional approach to making these prob- 
lems computationally tractable is to replace the Liou- 
ville/Kolmogorov equation by a (small) set of equations 
(PDEs or ODEs) for a few, low order statistical moments 
of its solution. When taking this approach for nonlinear 
systems, one must make an approximation, a closure, for 
the dependence of higher order moments on lower order 
moments. Typically the form of the closure equation is 
based on expert knowledge, empirical data, and/or phys- 
ical insight. For example, in the superposition approxi- 
mation and its extensions 0] for dense liquids and plas- 
mas, both quantum or classical, one approximates third 
order moments as functions of second order moments. 
Moment closure methods of this type have been applied 
to a number of areas including fluid turbulence (see Q 
and references therein, and also the work of Chorin et. 
al.). Of course, as with any approximation strategy, the 
quality of the resulting reduced description depends on 
the approximations made - poor closures lead to poor 
answers/predictions. In addition to replacing the ensem- 
ble with a small set of equations for low order moments, 
these equations are typically easier to solve. They are 
deterministic and generally far less stiff than the original 
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equations. 

A less exploited variant of this approximation scheme 
is the probability density function (PDF) based moment- 
closure approach. For PDF moment closures one makes 
an ansatz for the system statistics guided by available in- 
formation (e.g., symmetries). One then uses this ansatz 
in conjunction with the original dynamical equations to 
derive moment equations. Such PDF-based closures have 
been developed for reacting scalars advected by turbu- 
lence , phase-ordering dynamics 11 ] and a variety of 
other systems. This approach to moment-closure is a 
close analogue of the Rayleigh-Ritz method frequently 
used in solving the quantum-mechanical Schroedinger 
equation, by exploiting an ansatz for the wave-function. 
For a formal development of this point of view, see 12]. 

One of the obstacles to applying moment closures is 
that often the closure equations are too complicated to 
write down explicitly, even with the availability of com- 
puter algebra / symbolic computation systems. This 
is especially true for large-scale, complex systems, e.g. 
global climate models. Because of their great complexity, 
even if one could in principle derive the closure equations 
analytically, this procedure would be extremely difficult 
and time-intensive. Moreover, each time a model is up- 
dated, as climate and ocean models regularly are, the 
closure equations would have to be rederived. In other 
cases it may simply be impossible to determine the clo- 
sure equations analytically. This is especially likely when 
PDF's are not Gaussian, which is the case for most use- 
ful closures. Monte Carlo or other numerical methods 
may be needed in order to evaluate integrals for the mo- 
ments 13]. In addition, there may be situations where 
neither analytic nor numerical/MC integration will yield 
the closure equations due to the black-box nature of the 
available numerical simulator such as a compiled numer- 
ical code with an inaccessible source. Clearly, a need ex- 
ists for a robust approach to the general closure protocol 
which circumvents analytical difficulties. 

We address that need here by combining PDF closures 
with equation free modeling 1J] 15| . The basic premise 
of the equation- free method is to use an ensemble of short 
bursts of simulation of the original dynamical system to 
estimate, on demand, the time-evolution of the the clo- 
sure equations that we may not explicitly have. The 
equation-free approach extends the applicability of sta- 
tistical closures beyond the rare cases where they can be 
expressed in closed form. This hybrid strategy may be 
faster than the brute-force solution of a large ensemble of 
realizations of the dynamical equations since the closure 
version is generally smoother than the original problem. 

This paper is organized as follows. In Section 2 we 
describe the general features of PDF-based moment clo- 
sures. In Section 3 we explain how to implement the 
equation-free approach with these closures. We then, in 
Section 4, apply these ideas for a specific dynamical sys- 
tem, the stochastic Ginzburg-Landau (GL) equations us- 



ing a particular PDF-based closure scheme, the entropy 
method of Eyink and Levermore (ljj . We conclude with 
a discussion of closure quality, computational issues, and 
the application of our approach to large-scale systems. 



PDF-BASED MOMENT CLOSURES 

We consider the very general class of dynamical sys- 
tems, including maps, formally represented by 



X = U(X(t),N(t),t) 



U t (X t ,N t ) 



(1) 



(2) 



where N(i) is a stochastic process with prescribed statis- 
tics. The stochastic component arises from unknown pa- 
rameters, random forcing, neglected degrees of freedom 
and/or random initial conditions. This class includes 
both deterministic and stochastic systems with discrete 
and/or continuous states. Queueing systems, molecular 
dynamics, and stochastic PDEs are just some of the many 
examples that fall into this category. 

For concreteness in this paper we restrict ourselves to 
a special case of equation ([2]), namely, situations where 
N(£) is a Markov process (Brownian motion, Poisson pro- 
cess, etc.) and — more specifically still — Ito stochastic dif- 
ferential equations of the form: 



dX = U(X, t)dt + V2S(X, t)dW(t). 



(3) 



The deterministic component of the state, X, is gov- 
erned by the continuously differentiable vector field, U : 
R N xl-t R N . For many problems of interest (e.g., cli- 
mate) U is a highly nonlinear function. The noise com- 
ponent is modeled by the standard mean 0, covariance 
matrix I Wiener process, W 6 R , possibly modulated 



by a state-dependent matrix S 



tNxN 



Equa- 



tion ^ encompasses a wide class of systems including 
deterministic (S = 0) ones. 

In many cases one is interested in knowing the low 
order statistics of equation ([3]), for example an instanta- 
neous mean value or possibly multi-point covariance of 
X. These statistics can be obtained by averaging over 
an ensemble of stochastic systems, solving equation ([3]). 
They can also be obtained via the forward Kolmogorov 
equation for the probability density function P(X, i): 



dtP = C*(t)P, 



(4) 



where P satisfies the conditions: -P(X, t) > 0, and 
J -P(X, t) d~K = 1, and where C* is the generator of the 
Markov process. In the case of equation ([3]) this operator 
takes the form 

£*(tty(X) = -Vx-(U(X,^(X))+V^ : (D(X,i)V(X)). 

(5) 
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The forward Kolmogorov equation then becomes a 
Fokker-Planck equation 



8 t P + V x • (UP) = V x : (DP) 



(6) 



where D(X,t) = S(X, i)S(X, t) T is the nonnegative- 
definite diffusion matrix arising from the noise term. Un- 
like the original dynamical equation ([3]), the forward Kol- 
morogov equation (FKE) is both linear and determinis- 
tic. Dealing with it, therefore, has apparent advantages 
over the original ensemble of stochastic systems simu- 
lations. The price to pay for these advantages is that 
the FKE lives in a typically high, potentially infinite- 
dimensional, space. When equation |3]) is a nonlinear 
PDE, numerical solution to the FKE is usually ruled out. 

For computational purposes, we would therefore like 
to reduce the FKE (if possible and useful) to a small 
system of ordinary differential equations. This reduc- 
tion should simplify the computation as much as pos- 
sible while retaining fidelity to the original dynamical 
processes. The reduction proceeds by taking moments of 
the FKE with respect to a vector- valued function £(X, t) 
from R w x R + -> R M . The £ selected should include 
the relevant variables in the system (slow modes, con- 
served quantities, etc.). The moments fi(t) of £(X, t) are 
defined by 



M(t)= / C(X,t)P(X,t)dX 
and give rise to 



where 



C(x,f) = fit€(x,t) + A*)€(x,t) 



(7) 



(8) 



(9) 



and C is the adjoint of C* or the backward Kolmogorov 
operator. The result (jSJ) can be obtained by averaging 
over an ensemble of realizations of the stochastic dynam- 
ics ©. In general, however, © is not a closed equa- 
tion for the moments, fi. One can close this equation by 
choosing a PDF, P(X, t, f/,), which itself is a function of 
the moments fi. 



A(t)=V(M,*)s £(X,t)P(X,t,n)dX. 



(10) 



Alternatively, one can select a family of probability den- 
sities P(X, t, a), specified by parameters a — a(fi,t) 
rather than directly by the moments fi. This is analogous 
to specifying the temperature in the canonical ensemble 
as opposed to the average energy. The equivalence of 
these approaches is guaranteed provided that the param- 
eters and moments can be determined uniquely from one 
another. The translation between the parameters and 
their corresponding moments can be carried out by one 



of several methods. In some cases one may require Monte 
Carlo evaluation of the resulting integrals. 

If the moments and/or parameters are selected 
judiciously, one hopes that the approximate PDF 
P(X, t, ac(fi)(t)) will be close to the exact solution of 
the Liouville/Kolmogorov equation (g]). The mapping- 
closure approach of Chen et al llOl and the Gaussian 
mapping method of Yeung et al. I 111 are based on this 
type of parametric PDF closure [l9(. In fact, perhaps 
the most familiar application of the parametric approach 
is the use of the Rayleigh-Ritz method in quantum me- 
chanical calculations. This is the essential approach of 
our paper. 



EQUATION-FREE COMPUTATION 

Although we now have obtained a closed moment equa- 
tion (equation I10|) . we still need to determine the dynam- 
ical vector field V. As explained above, this step can be a 
serious obstacle to the practical implementation of PDF- 
based moment-closure (PDFMC). A method to calculate 
V is desirable that (i) does not require a radical revision 
each time the underlying code or model changes, and (ii) 
is relatively insensitive to the complexity of the PDFMC. 
The equation-free approach of Kevrekidis and collabora- 
tors 



14j meets those requirements. It permits one to 



work with much more sophisticated, physically realistic 
closures. 

Equation-free computation is motivated by the simple 
observation that numerical computations involving the 
closure equations ultimately do not require closed for- 
mulae for the closure equations. Instead, one must only 
be able to sample an ensemble of system states X dis- 
tributed according to the closure ansatz P(X, t; a) and 
then evolve each of these via equation ([3]) for short inter- 
vals of time. Such sampling and subsequent dynamical 
evolution would be necessary to calculate the statistics 
of interest even when not using a closure strategy. It is 
sufficient to have a (possibly black-box) subroutine avail- 
able which, given a specific state variable X(i) as input, 
returns the value of the state X(i + St) after a short 
time St. The ensemble of systems, each of which satis- 
fies equation ([3]), is evolved over a time interval St. The 
moments/parameters (iora are determined at the be- 
ginning and end of this interval and the time derivative 
(jl is estimated from the results of these short ensemble 
runs. This "coarse timestepper" can be used to estimate 
locally the right hand side of the closure evolution equa- 
tions, namely V(/x,i). 

Coarse projective forward Euler (arguably the simplest 
of equation-free algorithms) which we will use below il- 
lustrates the approach succinctly: Starting from a set 
of coarse-grained initial conditions specified by moments 
fj,(t) we first (a) lift to a consistent fine scale descrip- 
tion, that is, sample the PDF ansatz P(X, t;at(t)) to 
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generate ensembles of initial conditions X for equation 
d3j) consistent with the set £t(i); (b) starting with these 
consistent initial conditions we evolve the fine scale de- 
scription for a (relatively short) time St; we subsequently 
restrict back to coarse observables by evaluating the mo- 
ments fi(t + St) as ensemble- averages and (d) use the re- 
sults to estimate locally the time derivative dfi/dt. This 
is precisely the right hand-side of the explicitly unavail- 
able closure, obtained not through a closed form formula, 
but rather through short, judicious computational exper- 
iments with the original fine scale dynamics/code. Given 
this local estimate of the coarse-grained observable time 
derivatives, we can now exploit the smoothness of their 
evolution in time (in the form of Taylor series) and take 
a single long projective forward Euler step: 



Hit + At) = fi(t) + At 



H(t + St) - fi(t) 



St 



(11) 



The procedure then repeats itself: lifting, fine scale evo- 
lution, restriction, estimation, and then (connecting with 
continuum traditional numerical analysis) a new for- 
ward Euler step. Beyond coarse projective forward Eu- 
ler, many other coarse initial-value solvers (e.g. coarse 
projective Adams-Bashforth, and even implicit coarse 
solvers) have been implemented; the stability and accu- 
racy study of such algorithms is progressing These 
developments allow us to construct a nonintrusive imple- 
mentation of PDF moment closures, nonintrusive in the 
sense that we compute with the closures without explic- 
itly obtaining them, but rather by intelligently chosen 
computational experiments with the original, fine-scale 
problem. 

There is, however, an obvious objection to the 
equation-free implementation of moment-closures. Using 
the same ingredients, one can clearly obtain an estimate 
of any statistics of interest (for example, the moment- 
averages fJ,{i)) without the need of making any closure 
assumptions whatsoever. This can be done by the much 
simpler method of direct ensemble averaging. That is, 
one can sample an ensemble of initial conditions X from 
any chosen distribution P (X) , evolve each of these real- 
izations according to the fine-scale dynamics of equation 
(|3|), and then evaluate any statistics of interest at time 
t by averaging over the ensemble of solutions X(i). It 
would seem that this direct ensemble approach is much 
more straightforward and accurate than the equation- free 
implementation of a moment-closure, which introduces 
additional statistical hypotheses. 

The response to this important objection is that the 
fine-scale dynamics fl3J) is often very stiff for the appli- 
cations considered, in which the system contains many- 
degrees-of-freedom interacting on a huge range of length- 
and time-scales. In contrast, the closure equation (fT0|) is 
much less stiff, because of statistical-averaging, and its 
solutions n{t) are much smoother in time (and space). 
Thus, to evolve an ensemble of solutions of the fine-scale 



dynamics ([3]) from an initial time to to a final time to + T 
would require 0(T/St) integration steps, where the time- 
step St is required to be very small by the intrinsic stiff- 
ness of the micro-dynamics. In the closure approach, the 
evolution of the moment equations (|10p from time to to 
time to + T requires only 0{T / At) integration steps, with 
(hopefully) At ^S> St. Each of these closure integration 
steps by an increment At requires in the equation-free 
approach just one (or just a few) fine-scale integration 
step by an increment St. Thus, there is an over-all savings 
by a (hopefully) large factor 0(At/St). This crude esti- 
mate is based on a single step coarse projective forward 
Euler algorithm; clearly, more sophisticated projective 
integration algorithms can be used. 

In all of them, however, the computational savings are 
predicated on the smoothness of the closure equations, 
and are governed by the ratio of the time that it takes 
to obtain a good local estimate of dfi/dt from full direct 
simulation to the time that we can (linearly or even poly- 
nomially) extrapolate fj,(i) in time. It is also worth not- 
ing that a variety of additional computational tasks, be- 
yond projective integration (e.g. accelerated fixed point 
computation) can be performed within the equation-free 
framework 

In the next section we show by a concrete example how 
significant computational economy can be achieved with 
statistical moment closures implemented in the equation- 
free framework. 



A NUMERICAL EXAMPLE 

We illustrate here the equation-free implementation 
of moment-closures for a canonical equation of phase- 
ordering kinetics [It]], the stochastic time-dependent 
Ginzburg-Landau (TDGL) equation in one spatial di- 
mension. This is written as 



d<j>(x, i) 
dt 



DA<j>{x,i)-V'{<j>{x,t))+r){x,t) (12) 



where 4>(x, t) represents a local order parameter, e.g. a 
magnetization. The noise has mean zero and covariancc 
(rj{x,i)rj(x',t')) = 2kT5{x-x')8{t-t'). The potential V 
shall be chosen as 

to represent a single quartic/quadratic well. This 
stochastic dynamics has an invariant measure which is 
formally of Hamiltonian form P*[(f>] oc exp(— H [4>]/kT) 
where 



HI 



:^iv<m*)i 2 



V(<f>(x))]dx. (13) 



The Gibbsian measure P» [4>] is approached at long times 
for any random distribution Po[<p] of initial states. 
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One of the simplest dynamical quantities of interest is 
the bulk magnetization </>(t) = (l/V) J 4>(x,t)dx, where 
V is the total volume. If the initial statistics are space- 
homogeneous, then the ensemble average fi(t) = (<f>(t)) is 
also given by fi(t) = (4>(x, t)) for any space point x. Equa- 
tion ([T2]) leads to a hierarchy of equations for statistical 
moments of 4>(x, t). For example, the first moment satis- 
fies the equation 



dt 



A(<j>(x,t)) - (4>{x,t)) - (4> 3 (x,t)). (14) 



The evolution of the mean total magnetization is thus 
a function of the mean cubic total magnetization. One 
could write a time evolution equation for (</> 3 ), but it 
would involve a higher order term (<fi 5 ), and so on. Each 
equation contains higher moments and therefore the hi- 
erarchy does not close. 

To close the equation for /i(t) we assume a parametric 
PDF of the form P[4>] a] oc exp(— H[4>] a]/kT) where 



4>(x) dx 



is a perturbation of the Hamiltonian <| 1 3|) by a term pro- 
portional to the moment variable £[<f>] — (l/V) J 4>(x) dx. 
This is a special case of a general "entropy-based" clo- 
sure prescription proposed by Eyink and Levermore [l6| . 
This closure scheme guarantees that a(t) — > at long 
times and therefore the PDF ansatz P[<f>; ct(t)} relaxes to 
the correct stationary distribution P* [<f>] of the stochastic 
process. The determination of the parameter a given the 
moment jj, is here accomplished by Legendre transform 



a = argmax Q [a/i — F(a) 



(15) 



where the "moment- generating function" F(a) = 
log(exp[a J (fi(x)dx]}* and (•)* denotes average with re- 
spect to the invariant measure P* [</>] . The numerical op- 
timization required for the Legendre transform is well- 
suited to gradient descent algorithms such as the conju- 
gate gradient method, since 

(d / da)[an — F(a)] = fi — fi(a), 

where /J,(a) = is the average of the moment-function 
in the PDF ansatz P[<p;a\. In simple cases, F(a) and 
fJ,(a) — F'(a) may be given by closed analytical expres- 
sions. If not, then both of these averages may be deter- 
mined together by Monte Carlo sampling techniques. 

In the numerical calculations below, we discretize equa- 
tion (TT2]) using a forward Euler-Maruyama stochastic in- 
tegrator and 3-point stencil for the Laplacian (other dis- 
cretizations are possible). 

<f>(x, t + St) = (f>(x, t) - 5t[<p(x, t) + (/) 3 (x, t)] + (16) 
DSt 



(Sx) 



„ ^^x + Sx, t) — 2<f)(x, t) + 4>(x — Sx, t)} + 



where N(x, t) are independent, identically distributed 
standard normal random variables for each space-time 
point (x,t). The invariant distribution of the stochastic 
dynamics space-discretized in this manner has a Gibbsian 
form oc exp(—Hs/kT) with discrete Hamiltonian 



D 



2Sx 



(17) 



(x,x f ) 



Y J Sx[U\x) + -cp\x)] 



v / 2kT(St/Sx)N(x, t) 



where (x, x') are nearest-neighbor pairs. The closure 
ansatz can be adopted in the consistently discretized form 
Ps[4>;a] oc exp(— Hg[4>; a]/kT) where 

Hs[4>;a] = Hs[4>] + a^^Sx (j)(x). 



In this numerical experiment, we integrate an N = 
1000 member ensemble of solutions of equation (|17[) . 
and measure the ensemble-averaged, global magnetiza- 
tion fi(t) = (0(t)) = (l/V)J2 x ((/>(x,t)) at each time- 
step. With this we compare the results of the entropy- 
based closure simulation implemented by the equation- 
free framework using also an ensemble with N = 1000 
samples. In this concrete example, the projective inte- 
gration scheme works as follows: Suppose we are given 
the parameter a(t) at time t. The mean fj,(t) is first cal- 
culated from the parametric ensemble at time t by Monte 
Carlo sampling. Next all N samples are integrated over 
a short time-step St to create a time-advanced ensemble. 
From this ensemble fi(t + St) is calculated, which yields 
an estimate of the local time derivative. 

flapp(t) = [n(t + St)~ fi{t)]/6t. 

A large, projective Euler time-step of the moment aver- 
age is then taken via 

H(t + At) = n(t) + Atfi app (t). 

The parameter is finally updated by using the Legendre 
transform inversion to obtain a(t + At) from the known 
value n(t + At). The cycle may now be repeated to in- 
tegrate the closure equations by successive time-steps of 
length At. 

A critical issue in general application of projective inte- 
gration is the criterion to determine the projective time- 
step At. For stiff problems with time-scale separation, 
the projective time step for stability purposes is of the 
order of (1 /fastest "slow group" eigenvalues), while the 
"preparatory" simulation time is of the order of (1/slow- 
est "fast group" eigenvalue). Variants of the approach 
have been developed for problems with several gaps in 
their spectrum [18[ . Accuracy considerations in real-time 
projective step selection can, in principle, be dealt with 
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in the traditional way for integrators with adaptive step- 
size selection and error control: through on-line a poste- 
riori error estimates. An additional "twist" arises from 
the error inherent in the estimation of the (unavailable) 
reduced time derivatives from the ensemble simulations; 
issues of variance reduction and even on-line hypothe- 
sis testing (are the data consistent with a local linear 
model?) must be considered. These are important re- 
search issues that are currently explored by several re- 
search groups. Nevertheless, the main factor in computa- 
tional savings comes from the effective smoothness of the 
unavailable closed equation: the separation of time scales 
between the low-order statistics we follow and the higher 
order statistics whose effect we model (and, eventually, 
the time scales of the direct simulation of the original 
model) . 

Figure 1 is a plot comparing Projective Integration 
with Entropy Closure and direct Ensemble Integration 
with equation (p~2|) for diffusion constant D — 1000.0 We 
have selected both the "fine-scale" integration step St and 
the "coarse-scale" projective integration step At to be as 
large as possible, consistent with stability and accuracy. 
Thus, only steps small enough to avoid numerical blow- 
ups were considered. Then, values were selected both 
for St and for At so that the numerical integrations with 
those time-steps differed by at most a few percent from 
fully converged integrations with very small steps. In this 
manner, the time step required for the Euler-Maruyama 
integration of (fT2")l was determined to be St — 0.0004. On 
the other hand, for projective integration of the closure 
equation a time step At = 0.01 could be taken. This 
indicates a gain in time step by a factor of 25, which is 
also roughly the speed-up in the algorithm or savings in 
CPU time. The present example is not as stiff as equa- 
tions that appear in more realistic applications, with a 
very broad range of length- and time-scales, where even 
greater computational economies might be expected. 

In general, the moment-closure results need not agree 
so well with those of the direct ensemble approach, even 
when both are converged. In the example presented here, 
there is good agreement because the closure effectively 
captures the one-point PDF (see Fig. 2). This one-point 
PDF is the only statistical quantity that enters into Equa- 
tion (|14|) as long as the statistics are homogeneous and 
the Laplacian term vanishes. 

CONCLUSIONS 

In this paper, we have described how one can combine 
recently developed equation-free methods with statisti- 
cal moment closures to model nonlinear problems. With 
this method we can numerically integrate complex non- 
linear systems, for which closure equations may not be 
available in closed form. In the example presented here 
the specific entropy-based closure we selected has an H- 



80 




10 1 1 1 1 1 1 1 1 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 



FIG. 1: Mean total field as a function of time. Line (sym- 
bols): traditional (coarse projective) integration, respectively. 
See the text for a description of the stepsize selection. 
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FIG. 2: Comparison of the time dependent PDF's of the local 
field 4>(x, t) for the exact solution (blue) and for the projective 
integration / closure solution (red). 



theorem which guarantees relaxation to the equilibrium 
state of the original dissipative dynamics. However, we 
stress that the general approach outlined above can be 
used with a variety of closure methods. 

The equation-free method has the potential to enhance 
the flexibility, power, and applications set of the statis- 
tical moment closure approach. Since little or no an- 
alytic work is required, the sophistication of statistical 
moment closures can greatly enhanced beyond Gaussian 
PDF ansatze. The "practical usefulness" criterion for 
parametric PDF models that they permit analytical cal- 
culations is replaced by the criterion that they can be 
efficiently sampled. We believe that this approach can 
significantly increase the usefulness of closure methods. 

In order to model systems like global climate, oceans, 
and reaction diffusion processes in systems biology, one 
will have to construct more complex closures. These 
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will likely include higher order moments, correlation 
functions of the relevant variables, highly non-Gaussian 
statistics, etc. As the closures become more complex, 
the lifting step will require more efficient sampling ap- 
proaches. One will likely have to use nonlocal, acceler- 
ated sampling methods. One will also likely employ the 
latest in adaptive time and adaptive mesh methods to 
optimize performance for large-scale problems. 
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